■  X 

AD-A114  535  WISCONSIN  UNIV-MADtSON  MATHEMATICS  RESEARCH  CENTER  F/Q 

mathematical  M0DELIN6  OF  TWO-PhASE  FLOW.(U) 

MAR  82  D  A  DREW  DAA629-80-C-0041 

UNCLASSIFIED  MRC-TSR-2343  ML 


rnpy 


Mathematics  Research  Center 
University  of  Wisconsin— Madison 

610  Walnut  Street 
Madison,  Wisconsin  53706 

March  1982 


(Received  February  11,  1982) 


IjSponspred  by 


Approved  for  public  release 
Distribution  unlimited 


OTIC 


ELECTE 
MAY  1  8  1982 


S.  Army  Research  Office 
P.  0.  Box  12211 


ttesearch  Triangle  Park 
f^orth  Carolina  27709 


E 

82  G5  13  028 


UNIVERSITY  OP  WISCONSIN  -  MADISON 
MATHEMATICS  RESEARCH  CENTER 


MATHEMATICAL  MODELING  OF  TWO-PHASE  FLOW 
D.  A.  Drew 

Technical  Stanmary  Report  #  2343 

March  1982 
ABSTRACT 

A  continuum  mechanics  approach  to  two-phase  flow  is  reviewed.  An 
averaging  procedure  is  discussed  and  applied  to  the  exact  equations  of 
motion.  Constitutive  equations  are  supplied  and  discussed  for  the  stresses, 
pressure  differences  and  the  interfacial  force.  The  nature  of  the  resulting 
equations  is  studied. 
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SIGNIFICANCE  AND  EXPLANATION 


Equations  describing  the  motions  of  two  materials,  one  dispersed 
throughout  the  other,  are  derived.  There  is  some  controversy  over  the 
ability  of  certain  simplified  forms  of  these  equations  to  predict  meaning¬ 
ful  flows.  Consequently,  evidence  is  presented  to  verify  the  forms  of 
various  terms  appearing  in  the  equations.  The  consequences  of  the  assumed 
forms  are  discussed. 
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tlATHEMATICAL  MODELING  OF  TWO-PHASE  FLOW 
D.  A.  Drew 

INTRODUCTION 

Dispersed  two-phase  flows  occur  in  many  natural  and  technological 
situations.  For  example,  dust  in  air  and  sediment  in  water  contribute 
to  ertsion  and  silting,  and  can  cause  problems  for  machinery  such  as 
helicopters  and  power  plants  operating  in  such  an  environment.  Also, 
many  energy  conversion  and  chemical  processes  involve  two-phase  flows. 
Boiling  a  fluid  such  as  water  or  sodium  and  extracting  the  heat  by 
condensation  at  a  different  location  provides  a  practical  energy  flow 
process.  Many  kinds  of  chemical  reactants  and  catalysts  are  mixed  in  a 
dispersion  of  fine  particles  to  expose  as  much  interfacial  area  as 
possible. 

In  such  a  large  class  of  problems,  many  diverse  mechanisms  are 
important  in  different  situations.  The  models  used  for  these  different 
situations  have  many  common  features,  such  as  inter facial  drag. 

This  paper  examines  the  common  features  of  dispersed  two  phase 
flows  from  a  continuum  mechanical  approach.  Since  it  is  not  universally 
accepted  that  such  an  approach  is  valid,  we  shall  discuss  some 
philosophical  reasons  for  taking  such  an  approach.  The  approach  is 
based  on  the  view  that  it  is  sufficient  to  describe  each  material  as  a 
continuum,  occupying  the  same  region  in  space.  This  new  "material" 
consists  of  two  interacting  materials  (called  phases,  even  though  they 
often  are  not  different  phases  of  the  same  material).  The  two  phase 
material  is  often  called  the  mixture.  In  analogy  with  continuum 
mechanics,  we  shall  have  to  specify  how  the  mixture  interacts  with 
itself.  In  ordinary  continuum  mechanics  (ignoring  thermodynamic 
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considerations)  that  requires  a  constitutive  equation  expressing  the 
stress  as  a  function  of  various  fields.  In  two-phase  flow  mecnanics,  it 
requires  specification  of  stresses  for  each  phase,  plus  a  relation  for 
the  interaction  of  the  two  materials. 

Researchers  who  do  not  subscribe  to  this  view  treat  the  mechanics 
of  the  two  materials,  plus  the  dynamics  of  the  interface  as  fundamental, 
and  use  derived  results  and/or  measurements  to  gain  the  understanding 
needed  for  their  particular  application.  It  is  instructive  to  consider 
the  approaches  to  gas  dynamics.  Most  scientists  believe  that  a  gas  is  a 
collection  of  many  molecules  which  move,  vibrate,  and  interact  in  a 
complex,  but  describable  way.  Indeed,  with  the  aid  of  large  computers, 
molecular  dynamics  has  made  great  strides  in  understanding  such 
phenomena  as  shocks  and  phase  transition,  among  others.  In  spite  of  the 
knowledge  of  the  "correctness"  of  this  model,  many  scientists  and 
engineers  use  a  continuum  model  for  gas  dynamics.  Indeed,  anyone 
wishing  to  describe,  for  example,  the  flow  around  an  airfoil,  would  be 
hard  pressed  to  find  a  computer  big  enough  to  do  it  as  a  problem  in 
molecular  dynamics.  As  a  problem  in  continuum  gas  dynamics,  it  is  still 
a  large  problem;  however,  it  is  done  quite  routinely  numerically, 
including  shocks.  In  addition,  certain  solutions  to  the  continuum 
equations  can  be  obtained  analytically.  While  these  are  not  always  of 
direct  technical  interest,  they  often  suggest  phenomena  or  techniques 
which  do  have  direct  bearing  on  problems  of  interest.  Experiments  in 
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gas  dynamics  can  be  useful.  Without  a  model  for  comparison  or  scaling, 
however,  one  is  limited  to  understanding  only  the  particular  geometry 
and  scale  of  the  experiment.  The  ability  to  infer  is  lost. 

While  these  points  seem  obvious  for  gas  dynamics,  they  nevertheless 
should  be  discussed  and  understood  for  the  analogous  situation  of  two- 
phase  flow  mechanics.  For  a  particular  flow  with  ''nly  a  few  particles 
which  interact  only  slightly,  it  is  best  (and  perhaps  even  necessary)  to 
describe  them  “molecular ly " ,  that  is,  by  predicting  the  trajectories  of 
each  one.  If  many  particles  are  involved,  it  is  better  to  use  a 
continuum  description.  As  in  gas  dynamics,  both  descriptions  have  their 
place.  Furthermore,  it  is  often  instructive  to  try  to  ascertain  what 
one  model  implies  about  the  other,  as  in  using  the  Boltzmann  description 
of  a  gas  to  get  continuum  properties. 

The  analogy  is  apt  for  a  continuum  description  of  two-phase 
flows.  This  paper  will  review  the  connection  with  the  exact,  or 
microscopic  description  through  the  application  of  an  averaging  process 
to  the  continuum  mechanical  equations  describing  the  exact  motion  of 
each  material  at  each  point.  If  the  exact  flows  were  known,  the 
averaged  equations  would  be  completely  determined,  and  therefore 
unnecessary,  since  desired  averaged  information  (such  as  the  average 
concentration  of  particles)  could  be  determined  without  using  the 
averaged  equations.  The  solutions  are  not  known,  nor  is  it  necessary  to 
the  continuum  approach  that  they  be  known.  The  resulting  averaged 
e juations  are  assumed  to  describe  a  material  (the  mixture)  for  which  the 
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interactions  must  be  specified.  This  specification  is  then  done 
according  to  certain  rules  which  are  reasonable  and  not  too  limiting. 

We  shall  discuss  these  rules  and  their  implications. 

The  forms  of  the  resulting  equations  are  determined  by  the  choice 
of  a  list  of  variables  which  are  assumed  to  influence  the  interactions. 
The  resulting  equations  have  several  unknown  coefficients.  These 
coefficients  are  assumed  to  be  determinable  by  experiments.  For  the 
model  we  propose,  we  examine  the  experimental  data  and  their 
implications  on  the  coefficients.  The  result  is,  in  essence,  a 
recommendation  for  a  model  which  has  many  known  features  of  two-phase 
flow  dynamics.  An  enlightened  investigator  can  use  the  model  to  make 
predictions  in  a  situation  which  falls  within  the  range  of  assumptions 
made.  The  model  can  also  be  a  starting  point  to  obtain  generalizations 
(such  as  the  inclusion  of  electromagnetic  effects).  As  with  all  models, 
it  should  never  be  used  blindly,  but  with  caution  and  careful 
examination  of  the  results  and  implications. 

Historically,  geophysical  flows  involving  sediment  and  clouds  were 
among  the  first  two-phase  flows  observed  from  a  scientific  point  of 
view.  Sediment  meant  erosion  and  loss  of  property,  or  deposition  and 
loss  of  navigability.  Clouds  signalled  possible  rainfall.  The  amounts 
could  mean  drought,  bountiful  harvests  or  devastating  floods.  An 
analytical  description  awaited  the  development  of  the  calculus,  fluid 
mechanics  and  partial  differential  equation  theory.  The  industrial 
revolution  spurred  the  need  for  understanding  of  all  of  basic  science, 
including  two  phase  flows,  although  the  correlation  between  progress  and 


the  acquisition  of  fundamental  knowledge  is  not  exact.  For  example, 
automobiles  have  worked  reasonably  well  without  a  detailed  knowledge  of 
the  flow  and  evaporation  of  fuel  droplets;  on  the  other  hand,  efficient 
design  of  fluidized  beds  have  been  reliant  on  the  understanding  and 
stabilizing  effects  obtained  from  analytical  models. 

Early  work  on  the  form  taken  by  beds  of  particles  subject  to  forces 
due  to  flowing  fluid  (DuBuat  1786,  Helmholtz  1888,  Blasius  1912  and 
Exner  1920)  lead  to  interesting  fundamental  results  such  as  the 
stability  of  the  interface  between  two  fluids  and  a  practical 
understanding  of  the  macroscopic  properties  of  sedimentation  and  bed 
form  evolution.  A  more  microscopic  look  (Bagnold  1941)  dealt  with  the 
mechanics  of  the  interactions  between  the  particles  and  the  fluid,  and 
the  particles  and  the  bed.  A  recent  perspective  is  given  by  Engelund 
and  Fredsi^e  (1982). 

Developments  in  cloud  physics  occured  in  physical  chemistry  of 
nucleation  and  formation.  The  need  for  detailed  mechanical 
considerations  was  low  due  to  the  fact  that  the  vapor,  nucleation  sites 
and  the  formed  droplets  all  flow  with  the  surrounding  air  (until 
raindrops  form). 

Porous  medium  theory,  a  two  phase  flow  where  the  solid  phase  does 
not  move  appreciably,  developed  with  some  different  concerns  (Darcy 
1856),  but  Biot  (1955)  was  instrumental  in  putting  the  empirical 


knowledge  on  a  sound  continuum  basis.  The  desire  to  extract 


hydrocarbons  from  deep  inside  the  earth's  crust  has  given  a  great 
impetus  to  the  study  of  flow  in  porous  media.  See  Scheidegger 's  (1974) 
book  for  a  discussion  of  porous  bed  modeling  concerns. 

The  next  important  milestone  in  the  development  of  two-phase  flow 
continuum  mechanics  theory  was  the  development  and  use  of  equations  of 
dusty  gases  in  the  19S0's.  Models  for  dusty  gases  are  summarized  by 
Marble  (1970). 

Chemical  processing  in  fluidized  beds  gave  an  urgency  to  the 
development  of  the  theory  of  particle-fluid  systems.  Early  theoretical 
papers  include  those  of  Jackson  (1963),  Murray  (1965a,b),  and  Anderson 
and  Jackson  (1967,  1968).  The  flow  regime  involved  in  fluidization  is 
one  of  the  most  difficult  for  particle-fluid  flows.  The  particle 
concentrations  are  high,  the  dispersed  phase  is  relatively  dense,  the 
dispersed  phase  undergoes  a  random  micromotion,  and  often  is  is 
desirable  to  have  chemical  reactions  occur  in  the  flow. 

In  the  early  1960's  the  emergence  of  commercial  atomic  energy 
spurred  the  study  of  flows  of  steam  and  water.  The  work  of  Zuber  (1964) 
was  a  pioneering  landmark.  Fluid-fluid  flows  have  a  difficulty  not 
encountered  in  particle-fluid  flows,  namely  that  the  shape  of  the 
dispersed-phase  can  change,  resulting  in  changing  interfacial  area  and 
consequent  Interactions  between  the  fluids.  In  spite  of  much  progress 
(Lahey  &  Moody  1977),  two  phase  flow  studies  in  nuclear  reactors  are 


still  a  concern 


The  parallel  development  of  mixture  theory  had  some  small  influence 
on  the  progress  of  two-phase  flows.  Indeed,  the  concept  of 
interpenetrating  continuum  is  natural  in  mixtures  where  the  dispersion 
occurs  on  the  molecular  level.  The  development  of  the  ideas  of 
diffusion  by  Pick  (1855),  the  thermodynamic  concepts  in  mixtures  by 
Duhem  (1893),  Meixner  (1943),  Prigogine  and  Mazur  (1951)  lead  naturally 
to  the  theory  of  mixtures  expounded  by  Truesdell  and  Toupin  (1960).  The 
theory  of  mixtures,  as  applied  to  the  specific  mixtures  where  the  two 
constituents  remain  unmixed,  is  the  basis  of  the  present  work,  and 
(whether  explicitly  recognized  or  not)  much  of  the  more  specific 
previous  work.  Kenyon  (1976)  applies  the  mixture  theory  to  multiphase 
flows.  See  also  the  review  by  Bedford  and  Drumheller  (1982). 

By  the  1960's  enough  of  the  common  features  of  the  diverse  two- 
phase  flows  were  evident.  Several  influential  books  of  a  general  nature 
appeared,  including  Fluid  Dynamics  of  Multiphase  Systems  by  Soo  (1967), 
One-dimensional  Two-Phase  Flow,  by  Wallis  (1969),  and  Flowing  Gas-Solid 
Suspensions,  by  Boothrcyd  (1971).  Wallis'  approach  was  strongly 
influenced  by  gas-liquid  flows,  and  dealt  with  the  basic  concept  of 
interpenetrability  by  considering  cross-sectionally  averaged  equations, 
and  introduced  constitutive  assumption  by  quoting  appropriate 
experiments.  Soo's  work  was  largely  based  on  particle-fluid  flows.  He 
assumed  interpenetrability  from  the  start,  and  included  forces  in  the 
particular  momentum  equation  known  from  experiments  or  inferred  from 
calculations.  Boothroyd  was  interested  in  particle-gas  flows,  and 
contributed  ideas  on  turbulence  and  drag  reduction. 


A  different,  but  somewhat  related  approach  to  the  problem  of  a 
mechanical  description  of  two-phase  flows  uses  the  single-particle  or 
several  interacting  particle  flow  fields,  along  with  an  averaging 
approach  to  yield  rheological  or  transport  properties.  The  celebrated 
result  in  this  area  is  the  effective  viscosity  result  due  to  Einstein 
(1906),  which  shows  that  in  slow  flow  the  mixture  behaves  li)ce  a  fluid 
with  viscosity  increased  by  a  factor  of  1  +  ^  where  a  is  the 
volumetric  concentration  of  particles.  An  approach  outlined  by  Brinkman 
(1947)  has  been  influential.  The  philosophical  roots  of  this  approach 
are  elegantly  discussed  in  the  book  Low  Reynolds  Number  Hydrodynamics, 
by  Happel  and  Brenner  (1965).  The  idea  is  to  use  solutions  of  the  flow 
equations  in  special  cases  (such  as  Stoke 's  flow  around  an  array  of 
spheres)  to  infer  information  about  the  flow  (such  as  the  total  force  on 
the  particulate  phase).  Several  aspects  of  this  work  have  been 
discussed  in  this  Review;  see  the  papers  by  Brenner  (1970),  Batchelor 
(1974),  Herczynski  and  Pienkowska  (1980),  Leal  (1980)  and  Russel  (1981). 

Several  authors  (Bedford  &  Drumheller  1978,  Drumheller  &  Bedford 
1980)  have  pursured  a  variational  approach  to  two-phase  flows, 
generalizing  the  work  of  Biot  (1977).  The  variational  formulation 
starts  with  Hamilton's  principle,  given  by 

t,  t 

6  /  (T-V)dt  +  /  6wdt  =  0 
2  *^1 

where  6  represents  the  variation  over  an  appropriate  space  of 
functions,  T  and  V  are  the  kinetic  and  potential  energies,  6w  is 
the  virtual  work  and  t^  and  t2  are  two  arbitrary  times.  A 


variational  formulation  has  an  advantage  in  that  if  it  is  desired  to 
include  a  certain  effect  (for  example,  the  virtual  work),  that  effect 
would  be  included  consistently  in  the  mass,  momentum,  and  kinetic  energy 
equations.  A  more  concrete  and  immediate  advantage  is  in  formulating 
numerical  techniques,  and  specifically,  finite  element  techniques,  where 
variational  formultion  eases  the  translation  of  the  partial  differential 
equations  into  discrete  equations. 

In  order  to  use  a  variational  formulation,  it  is  necessary  to 
define  the  variation  <5,  the  kinetic  and  potential  energies  T  and 
V  and  the  virtual  work  6w.  Constraints  must  be  included  with  Lagrange 
multipliers.  Bedford  &  Drumheller  (1982)  note  that  Hamilton's  principle 
is  usually  formulated  with  a  "control  mass",  that  is,  for  material 
volumes.  The  problem  of  having  two  materials  moving  at  different 
velocities  is  dealt  with  by  assuming  a  control  volume  which  coincides 
with  a  rigid  surface  through  which  no  material  of  either  phase  passes. 

They  use  the  technique  to  show  how  the  effect  of  oscillations  in 
bubble  diameter  can  be  included  by  including  the  kinetic  energy  of  the 
liquid  due  to  a  change  in  bubble  diameter.  They  similarly  include 
virtual  mass  by  arguing  that  relative  accelerations  increase  the  kinetic 
energy.  A  conceptual  difficulty  with  the  approach  is  that  the  user  must 
decide  which  fluctuations  contribute  to  the  total  kinetic  energy,  and 
which  contribute  to  virtual  work.  For  example,  viscous  drag  is  due  to 
"fluctuations"  of  the  fluid  velocity  near  a  particle, but  is  included  as 
a  virtual  work  term.  A  rule  of  thumb  might  be  to  consider  whether  the 
energy  is  recoverable  or  not.  Virtual  work  is  associated  with 
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unrecoverable  energy,  kinetic  and  potential  energy  with  recoverable 
energy.  Variational  formulations  (in  general)  do  not  mention  the  effect 
of  Reynolds  stresses,  which  are  one  manifestation  of  fluctuations,  we 
conclude  that  while  variational  formulations  are  useful,  it  is  not 
always  straightforward  to  formulate  them. 
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EQUATIONS  OF  MOTION 

We  start  by  assuming  that  each  material  involved  can  be  described 
as  a  continuum,  governed  by  the  partial  differential  equations  of 
continuum  mechanics.  The  materials  are  separated  by  an  interface,  which 
we  assume  to  be  a  surface.  At  the  interface,  jump  conditions  express 
the  conditions  of  conservation  of  mass  and  momentum. 

The  equations  of  motion  for  each  phase  are  (Truesdell  and  Toupin 

1960) 

( 1 )  conservation  of  mass 

1^  +  7*pv  =0  (1) 

(2)  conservation  of  linear  momentum 

+  V*pw  =  7*t  +  pf  (2) 

valid  in  the  interior  of  each  phase.  Here  P  denotes  the  density,  v 
the  velocity,  T  the  stress  tensor,  and  £  the  body  force  density. 
Conservation  of  angular  momentum  becomes  T  =  T^,  where  t  denotes  the 
transpose.  At  the  interface,  the  jump  conditions  are 

( 1 )  jump  condition  for  mass 

( [P(v  -  v^)*n]  ]  =  0  (3) 

(2)  jump  condition  for  momentum 

[[Pv(v  -  v.)»n  -  T*n]]  =  axn  .  (4) 

1 

Here  [  [  ]  )  denotes  the  jump  across  the  interface,  is  that 

velocity  of  the  interface,  o  is  the  surface  tension,  assumed  to  be  a 
constant,  <  is  the  mean  curvature  of  the  interface,  and  n  is  the 
unit  normal  (Aris  1962).  We  shall  assume  that  n  points  out  of  phase 
k,  and  that  the  jump  between  f  in  phase  k  and  f  in  phase  £  is 
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i  k 

defined  by  [[fll  =  f  -  f  ,  where  a  superscript  k  denotes  the 
limiting  value  from  the  phase  k  side.  As  a  sign  convention  for  the 
curvature,  we  assume  that  <  is  positive  (concave)  toward  -n.  The 
mass  of  the  interface  has  been  neglected,  and  the  surface  stresses  have 
been  assumed  to  be  in  the  form  of  classic  surface  tension. 

We  do  not  discuss  any  thermodynamic  relations  in  this  paper. 
Thermodynamic  considerations  are  important  for  many  multiphase  flows. 

Our  discussion  focuses  on  the  mechanics  of  two-phase  flows,  hence  we 
elect  to  forego  a  discussion  of  thermodynamics  for  the  sake  of 
simplicity. 

Constitutive  equations  must  be  supplied  to  describe  the  behavior  of 
each  material  incolved.  For  example,  if  one  material  is  an 
incompressible  liquid,  then  specifying  the  value  of  p,  and  assuming 
T  =  -pi  +  p (Vv  +  (7v)^)  determines  the  nature  of  the  behavior  of  the 
fluid  in  that  phase.  Similar  considerations  are  possible  for  solid 
particles  or  a  gas.  The  resulting  differential  equations,  along  with 
the  jump  conditions,  provide  a  fundamental  description  of  the  detailed 
or  exact  flow. 

Usually,  however,  the  details  of  the  flow  are  not  required.  For 
most  purposes  of  equipment  or  process  design,  averaged,  or  macroscopic 
flow  information  is  sufficient.  Fluctuations,  or  details  in  the  flow 
must  be  resolved  only  to  the  extent  that  they  effect  the  mean  flow  (like 
the  Reynolds  stresses  effect  the  mean  flow  in  a  turbulent  flow). 
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In  order  to  obtain  equations  which  do  not  contain  the  details  of 
the  flow,  it  has  become  customary  to  apply  some  sort  of  averaging 
process.  It  is  not  essential  to  do  so;  indeed,  some  researchers  prefer 
to  postulate  macroscopic  equations  without  reference  to  any  microscopic 
equations.  Certainly,  the  necessary  terms  in  the  macroscopic  equations 
can  be  deduced  without  using  an  averaging  process.  One  advantage  of  a 
postulational  approach  is  obvious  -  not  having  to  deal  with  the  worries 
of  the  averaging  process.  The  advantages  of  averaging  are  less  obvious. 
First,  the  various  terms  appearing  in  the  macroscopic  equations  are 
shown  to  arise  from  appropriate  microscopic  considerations.  For 
example,  stress  terms  arise  from  microscopic  stresses  (pressure,  for 
example)  and  also  from  velocity  fluctuations  (Reynolds  stresses). 
Knowledge  of  this  fact  does  give  a  modicum  of  insight  into  the 
formulation  of  constitutive  equations.  If  a  term  appeared  in  the 
averaging  which  was  not  expected  in  the  postulational  approach,  then  it 
would  be  obligatory  to  include  it  in  the  postulated  model,  or  else 
explain  why  it  is  superfluous. 

An  additional  advantage  of  averaging  is  that  the  resulting 
macroscopic  variables  are  related  to  microscopic  variables.  If  the 
microscopic  problem  can  be  solved  for  some  special  situation,  the 
solution  can  then  be  used  to  get  values  of  the  macroscopic  variables. 
This  cannot  replace  the  need  to  postulate  constitutive  equations,  but  is 
can  give  insights  into  the  types  of  terms  expected  to  be  important  in 
tile  constitutive  relations.  The  philosophy  taken  here  is  that  whatever 


-14 


information  can  be  gained  from  averaging  is  worthwhile/  and  therefore  we 
present  a  generic  averaging  method/  and  its  results.  If  the  reader 
fee..3  that  averaging  is  unnecessary,  he  can  skip  this  section,  and 
assume  the  equations  (40)  -  (43)  have  been  postulated,  along  with  the 
interpretations  given  at  the  end  of  this  Section. 

Averaging  the  equations  of  motion  is  suggested  by  the  averaging 
approach  to  turbulence  (see  Hinze  1959).  Time  and  space  averages  were 
the  first  to  appear  (Frankl  1953,  Teletov  1958).  These  averages  have 
been  refined  by  weighting,  by  multiple  application,  and  by  judicious 
choice  of  averaging  region.  The  highlights  can  be  found  in  the  work  of 
Anderson  and  Jackson  (1967),  Vernier  and  Delhaye  (1968),  Drew  (1971), 
Whitaker  (1973),  Ishii  (1975),  Nigmatulin  (1979),  and  Gough  (1980). 
Statistical  averages  are  most  convenient  for  the  rheology  work;  see 
Batchelor  (1974)  for  a  summary.  The  application  of  statistical  averages 
to  the  equations  of  motion  is  straightforward;  the  paper  of  Buyevich  and 
Shchelchkova  (1978)  summarizes  the  approach  nicely. 

Let  <  >  denote  an  averaging  process  so  that  if  f(x,t)  is  an 
exact  microscopic  field,  then  <f>(x,t)  is  the  corresponding  averaged 
field.  We  shall  specify  shortly  the  requirements  which  an  average 
should  possess,  however,  for  now  we  merely  suggest  that  it  should  be 
smoothing  in  the  sense  that  no  details  appear  in  the  averaged 
variables.  Some  examples  of  commonly  used  averages  in  multiphase  flow 
are  the  time  average: 

<f>,(x,t)  -  ^  f(*,f)df  (5) 

where  T  is  an  averaging  time  scale;  the  space  average: 
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1  ^ 

<f>,(x,t, =^/ ’ :  / 


x^+  xL  x+“  L  x+“L 

^  ^  ^  ^  f (x’ ,t)dx^dx^dx^ 


(6) 


where  L  is  an  averaging  length  scale;  a  weighted  space  average 
<f>  (x,t)  =  ///  g(x-x' )f (x' ,t)dx' 


(7) 


where  ///  g(8)d8  =  1,  and  various  combinations  of  averages  and/or 
r3 

specific  types  of  weightings-  Also  mentioned  in  the  literature  are 
ensemble  averages: 

1  ^ 

<f>^(x,t)  =  -  E 
n=1 
or 


(8) 


<f>  (x,t)  =  /o  f (x,t;u»)dii(iij) 


(9) 


where  fj^(x,t)  or  f(x,t;(‘J)  denote  a  realization  of  the  quantity  f 
over  a  set  of  possible  "equivalent"  realizations  £1.  One  way  in  which 
randomness  might  be  introduced  in  a  particular  flow  situation  is  by 
allowing  the  particles  to  have  random  initial  positions. 

The  averaging  process  is  assumed  to  satisfy 


<f+g>  =  <f>  +  <g> 

(10) 

<<f>g>  =  <f><g> 

(11) 

<c>  =  c 

3f  3 

<3t>  =  ^ 

(12) 

3f  3 

'^3x.^  3x,  • 

1  1 

(13) 

The  first  three  of  these  relations  are  called  Reynolds  rules,  the 


fourth  is  called  Liebnitz'  rule,  and  the  fifth  is  called  Gauss'  rule. 


-16- 


Some  difficulty  is  encountered  when  trying  to  apply  the  average  to 
the  equations  of  motion  for  each  phase.  In  order  to  do  this,  we 
introduce  the  phase  function  Xj^(x,t)  which  is  defined  to  be 


X^(x,t)  = 


1  if  X  is  in  phase  k  at  time  t 


0  otherwise 


(14) 


We  shall  deal  with  X^  as  a  generalized  function,  in  particular  in 
regard  to  differentiating  it.  Recall  that  a  derivative  of  a  generalized 
function  can  be  defined  in  terms  of  a  set  of  "test  functions"  which 

are  "sufficiently  smooth"  and  have  compact  support.  Then  qr”  and 

T —  are  defined  by 
dx 


3X  g 

/  3  g—  (x,t)<t>(x,t)dxdt  =  -  /  ,  X^(x,t)  (x,t)dx(Jt  , 

R  XR 


R^R 


(15) 


3X 

/  j  g —  ( X,  t )  'K  X,  t )  dxdt  = 

R  XR  ’‘i 


-  /  3  Xj^(x,t)  -1^  (x,t)dxdt  .  (16) 

R  XR 


It  can  be  shown  that 

sT  *  ''1-%  ■  " 

in  the  sense  of  generalized  functions.  To  see  this,  consider 


3X 

^3  \*^xj^dxdt 

R  ^  R 


R^xr  ^ 


(17) 
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=  ■  /r« 

=  -  c  (ft’  4.  (t)* 

k 

=  /r  1_C0  =  °  ' 

since  if  has  compact  support  in  t.  Here  is  the  region 

occupied  by  phase  k  at  time  t,  and  we  assume  that  is  extended 

smoothly  through  phase  k  {in  order  that  the  second  line  makes  sense). 
If  f  is  smooth  except  at  S,  then  defined  via 

fVx^()idxdt  =  -  /  3  Xj^7(f()))d*dt 

R^xr  r'^xr 

'  <(-«  4^(t)  ^(^*)‘***’= 

=  r.  /s 

where  nj^  is  the  unit  normal  exterior  to  phase  k,  and  f  denotes 
the  limiting  value  of  f  on  the  phase-k  side  of  S, 

It  is  also  clear  that  is  zero,  except  at  the  interface. 

Equation  (19)  describes  the  behavior  of  Vx^  at  the  interface.  Note 
that  it  behaves  as  a  "delta-function",  picking  out  the  interface  S, 
and  has  the  direction  of  the  normal  interior  to  phase  k. 

This  motivates  writing 

=  \  Tn 


(20) 
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where  3x/3n  is  a  scalar  walued  generalized  function  which  has  the 
property  that 

/  -  (x,t)  ](>(*, t)d»dt  =  -/*/(;  ♦(x,t)dSdt  . 


(21 ) 


The  quantity  3x/3n  then  picks  out  the  interface  S.  We  write 

<|i>  =  s  (22) 

where  a  is  the  average  interfacial  area  per  unit  volume. 

Averaged  Equations 

In  order  to  derive  averaged  equations  for  the  motion  of  each  phase, 
we  multiply  the  equation  of  conservation  of  mass  valid  in  phase  k  (1) 
by  Xj^  and  average.  Noting  that 

3X. 


and 


we  have 


3p  3  k  3 

\  4  ■  TC  V  -  ^  ST-  ■  IT  ’St"  * 


7*pv  =  V»Xj^pv  -  pv*VXj^ 


3  )c 

<\f>>  +  V*<Xj^Pv>  =  <[p(v-v^)]^Vx^>  . 

Similar  considerations  for  the  momentum  equations  yield 


(23) 


(24) 


(25) 


<Xj^P^>  +  V«<Xj^pw>  =  V*<Xj^T>  +  <X^Pf> 


+  <[pv(v-v. 


i)  -  Tl-.VX^ 


>  . 


The  terras 


and 


<[p(v-ir  ) J^»VXj^>  = 


:(pv(v-v^)  -  T]^»VXj^>  = 


(26) 


(27) 


(28) 
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are  the  interfacial  source  terms.  As  noted,  Vx^  pic)cs  out  the 
interface,  and  causes  discontinuous  quantities  multiplying  it  to  be 
evalued  on  the  phase-k  side  of  the  interface. 

The  jump  conditions  are  derived  by  multiplying  equations  (2.3)  and 
(2.4)  by  3x/3n,  and  recognizing  that  We  obtain 


i  =  I  \  =  0 


k=1 


k=1 


I  <[pv(v-v.  )  -  T]^*Vx  >  =  0<icVx  > 

k=l  ^  k  1 


(29) 


(30) 


=  M 

m 

The  term  "m  is  the  contribution  to  the  total  force  or  the  mixture  due 
to  the  interface,  and  specifically  due  to  surface  tension. 

Applying  a  different  averaging  process  requires  a  different  set  of 
machinations  regarding  the  interfacial  source  terms  (Anderson  fi  Jackson 
1967,  Drew  1971,  Ishii  1975,  Delhaye  and  Achard  1979).  Almost  all  of 
the  derivations  for  specific  averaging  processes  seem  to  be  more 
complicated  than  the  above;  however,  the  trade-off  for  the  simple 
derivation  is  that  all  manipulations  now  involve  generalized  functions. 

The  volumetric  concentration  (or  volume  fraction,  or  relative 
residence  time)  of  phase  k  is  defined  by 

“k-<V  • 


We  note  that 


■3t'  " 


(32) 


and 


Vo,  =  <Vx  >  . 

k  k 


(33) 


I 
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There  are  two  types  of  averaged  variables  which  are  useful  in  two- 
phase  mechanics,  namely  the  phasic,  or  Xj^-weighted  average,  and  the 
mass-weighted  average.  Which  is  appropriate  is  suggested  by  the 
appearance  of  the  quantity  in  the  equation  of  motion.  The  phasic 
average  of  the  variable  is  defined  by 

and  the  mass  weighted  average  of  the  variable  is  defined  by 

=  <X.  P<Ji>/oi.  p.  .  (35) 

k  k  k  k 

It  is  convenient  to  write  the  stresses  in  terms  of  pressures 

plus  extra  stresses.  Thus, 

\  °  *  \  * 

It  is  expected  that  readers  familiar  with  fluid  dynamical  concepts  are 
familiar  with  the  concept  of  pressure  in  fluids;  in  this  case,  p^  can 
be  thought  of  as  the  average  of  the  microscopic  pressure.  If  one  of  the 
phases  consists  of  solid  particles,  the  concept  is  less  familiar.  In 
this  case,  the  microscopic  stress  (involving  small  elastic  deformations, 
for  example)  is  thought  of  being  made  up  of  a  spherical  part  (acting 
■equally  in  all  directions)  plus  an  extra  stress.  The  spherical  part, 
when  averaged,  yields  the  pressure  p^^  in  equation  (36). 

It  has  further  become  customary  to  separate  various  parts  of  the 
interfacial  momentum  transfer  term.  This  is  done  by  defining  the 
interfacial  velocity  of  the  k^^  phase  by 

T  v  .  =  <[pv(v-v.  )]*^*Vx  >  ,  (37) 

K  K  ^  X  X  K 

and  the  interfacial  pressure  on  the  k*"”  phase  by 


i 


Equation  (38)  is  the  dot  product  of  of  the  "standard"  definition 

(Ishii  1975)  of  the  interfacial  pressure.  The  standard  definition  uses 
three  equations  to  define  one  scalar  quantity,  and  cannot  be  a  generally 
valid  definition.  Here  the  remaining  part  of  the  contribution  of  the 
pressure  at  the  interface  is  lumped  with  the  viscous  stress  contribution 
at  the  interface,  and  is  treated  through  the  use  of  a  constitutive 
equation.  Thus,  we  write 

-  ■'kVi  - 

d  }c  )c 

where  M.  =  <(p“p,  .  )  -  V*^X  >  is  referred  to  as  the  interfacial 

K  pC  ^  X  jC  iC  iC 

force  density,  although  it  does  not  contain  the  effect  of  the  average 
force  on  the  interface  due  to  the  average  interfacial  pressure.  The 
term  -p^^  i^^k'  does  contain  the  force  due  to  the  average 

interfacial  pressure,  is  sometimes  referred  to  as  the  buoyant  force. 

The  reason  for  this  terminology  is,  of  course,  that  the  buoyant  force  on 
an  object  is  due  to  the  distribution  of  the  pressure  of  the  surrounding 
fluid  on  its  boundary. 

With  equations  (31)  and  (34)  -  (39),  the  equations  of  motion  (25) 


and  ( 26 )  become 


3a  P 

k  k  „  A  „ 

■3r“  ^  =  ^k 


^Vk\  „  -  A  A  ;  „  . 


^  Vk,i  ^  ^Pk,i  -  Pk^^“k 


(41  ) 
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The  jump  conditions  (29)  and  (30)  are 


I 


k=1 


0 


I 


)c=1 


^^k\.i  ^  Pk,i% 


M 

m 


(42) 


(43) 


Adequate  models  for  compressibility  and  phase  change  require 

consideration  of  thermodynamic  processes.  These  are  beyond  the  scope  of 

this  paper;  therefore  we  shall  restrict  our  attention  to  incompressible 

materials  where  no  phase  change  occurs.  Thus  we  assume  that 

p,  =  constant  (44) 

K 

and 

r,  =  0  .  (45) 

k 

In  order  to  simplify  the  notation,  we  shall  drop  all  symbols 
denoting  averaging. 
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CONSTITUTIVE  EQUATIONS 

In  order  to  have  a  useable  model  for  an  ordinary  single-phase 
material,  relations  must  be  given  which  specify  how  the  material 
interacts  with  itself.  These  relations  specify  how  the  material 
transmits  forces  from  one  part  of  the  body  to  another.  For  single-phase 
materials  a  relation  between  the  stress  and  other  field  variables  is 
usually  required.  For  two-phase  materials,  where  we  desire  to  track 
both  phases,  we  must  specify  the  stresses  for  both  phases  and  the  way  in 
which  the  materials  interact  across  the  interface.  Thus,  we  require 
constitutive  equations  for  the  stresses  (T  +  )  r  the  interfacial 

force  density  ,  and  the  pressure  differences  i'  consistent 

with  the  equations  of  motion  and  the  jump  conditions. 

The  fundamental  process  consists  of  proposing  forms  for  the 
necessary  terms  within  the  framework  of  the  principles  of  constitutive 
equations,  finding  solutions  of  the  resulting  equations,  and  verifying 
against  experiments.  The  process  can  be  iterative,  with  some  experiment 
indicating  that  a  more  involved  form  is  needed  for  some  constituted 
variable.  The  ideal  end  result  of  the  process  is  a  set  of  equations 
which  could  be  used  to  predict  the  behavior  of  the  two-phase  flow,  for 
example  with  a  computer  code.  With  the  equations  should  come  a  set  of 
conditions  for  the  validity  of  the  values  of  the  constants  and  other 
functions  used  in  the  constitutive  equations. 

The  stresses  +  o^,  the  interfacial  force  density  and  the 

pressure  differences  “  Pk  i  assumed  to  be  functions  of  a^, 

3oi  /3t,  ,  V  ,  ,  3v  /3t,  ...  ,  where  ...  represents  the  material 

K  K  K  K  K 


t 
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properties,  such  as  the  viscosities  and  densities  of  the  two  materials, 
and  other  geometric  parameters  such  as  the  average  particle  size,  or  the 
interfacial  area  density. 

For  concreteness,  we  shall  refer  to  phase  one  as  the  particulate, 
or  dispersed  phase  and  include  in  that  description  solid  particles, 
droplets,  or  bubbles.  Phase  two  is  then  the  continuous,  or  carrier 
phase,  and  can  be  liquid  or  gas.  We  shall  denote 

Cl  =  (46) 

so  that 

1-a  =  .  (47) 

It  is  evident  that  both  c»  and  1-a  need  not  be  included  as 
independent  variables  in  forming  constitutive  equations. 

Drew  and  Lahey  (1979)  consider  the  general  process  of  constructing 
constitutive  equations.  The  principle  of  consequence  is  that  of 
material  frame  indifference,  or  objectivity.  This  principle  requires 
that  constitutive  equations  be  appropriately  invariant  under  a  change  of 
reference  frame.  The  motivation  behind  this  assumption  is  that  we 
expect  the  way  in  which  a  two-phase  material  distributes  forces  to  be 
independent  of  the  coordinate  system  used  to  express  it.  Those  who 
doubt  this  principle  often  argue  that  the  momentum  equation  is  not  frame 
indifferent,  needing  Coriolis  forces  and  centrifugal  forces  to  correct 
it  in  non-inertial  frames.  This  argument  is  incorrect,  since  the 
appropriate  formulation  of  the  momentum  equation  deals  not  with 
"acceleration"  but  with  "acceleration  relative  to  an  inertial  frame". 
When  formulated  with  the  "acceleration  relative  to  an  inertial  frame". 
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the  momentum  equations  are  frame  indifferent  (Truesdell  1977).  This 
introduces  an  interesting  question:  Can  the  mixture  know  whether  or  not 
it  is  being  referred  to  a  inertial  frame?  That  is,  can  the  interfacial 
force  and  the  Reynolds  stresses  depend  on  the  reference  frame?  The 
principle  of  material  frame  indifference,  as  used  here,  claims  that  they 
cannot.  Others  claim  that  they  must  (Ryskin  and  Rallison  1981). 

The  principle  of  objectivity  is  now  examined  (Drew  s  Lahey  1979). 
The  approach  we  take  deals  with  coordinate  changes,  but  is  equivalent  to 

more  abstract  approaches  (Truesdell  1977).  Consider  a  coordinate  change 

_  * 

from  system  x  to  system  z  ,  given  by 

X  =  Q( t) ‘x  +  b(t)  ,  (48) 

where  Q(t)  is  an  orthonormal  tensor  and  b  is  a  vector.  A  scalar  is 
objective  if  its  value  is  the  same  in  both  the  starred  and  unstarred 
systems,  that  is,  if 

(j)  (x  ,t)  =  (P(x,t)  .  (49) 

A  vector  is  objective  if  it  transforms  coordinates  correctly,  that  is, 
if 

* 

▼  =  Q'v  .  (50) 

A  tensor  is  objective  if 

*  t 

T  =  Q*T*Q  .  (51  ) 

The  scalar  a  is  objective;  however  3a/3t  is  not.  It  is 

straightforward  to  show  that  D  a/Dt  =  3a/3t  +  v  • 7a  is  objective. 

k  k 

Let  us  next  consider  objective  vectors.  The  volume  fraction 
gradient  7a  is  objective.  If  we  differentiate  (48)  with  respect  to 
t  following  a  material  particle  of  phase  k,  we  see  that 


=  Q(t)*V|^  +  (Q»x  +  b) 


(52) 


Hence  velocities  are  not  objective.  This  is  obvious  physically,  since 
the  velocity  of  a  point  depends  on  the  coordinate  system.  If  we  take 
equation  (2.53),  for  k  =  1  and  2,  and  subtract,  we  have 

Thus,  the  relative  velocity  between  the  phases  is  objective. 

For  acclerations,  differentiating  (52)  following  a  material 
particle  of  phase  j  yields 


Thus, 

j  =  2 


D  .  V,  D  .  V,  .. 

i  k  i  k  *  * 

~  •  =  Q*  "rr~  +  C**  +  Q*v.  +  Q*v,  +  b  . 

Dt  *  Dt  *  “  3  *  k 


subtracting  the  value  for  j  =  1,  k  =  1  from  the  value  for 
k  =  1  yields  the  result  that 

D  V  D  V 

1  2  2  1 


t2  Dt 


Dt 


(54) 


(55) 


is  objective. 

The  list  of  tensors  which  we  consider  is  motivated  by  the  desire  to 
model  two-phase  flows.  Thus,  velocity  gradients  (expressing  rate  of 
deformation)  are  included,  but  displacement  gradients  (expressing 
deformation)  are  not.  It  is  well  known  (Truesdell  and  Toupin  1960)  that 

“k.b  - 1  <%  * 

are  objective.  We  further  note  that  ~  ^2  ^  objective 

tensor. 

Thus,  the  problem  of  forming  constitutive  equations  reduces  to 
expressing 


+  o 


Pv  -  Pv 


(57) 
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in  terms  of 


a,  Dj^o/Dt,  Va,  -  v^,  /Dt  -  v^/Dt, 


°1,b'  -  ^2’  • 


(58) 


We  shall  assume  that  the  materials  are  isotropic.  This  means  that 
no  direction  is  preferred  in  either  material  or  in  the  way  they 
interact.  Note  that  a  particular  flow  may  occur  in  such  a  way  that  a 
preferred  direction  (down,  for  example)  might  emerge.  In  that  case,  the 
flow  is  anisotropic,  but  the  materials  and  the  mixture  are  isotropic. 

If  f  is  a  scalar,  then  f  can  depend  on  the  scalars  and  the 
scalar  invariants  formed  from  the  vectors  and  tensors  in  the  list 
(58).  The  length  of  a  vector  is  a  scalar  invariant,  as  are  the  trace, 
determinant  and  Euclidean  norm  of  a  tensor. 

If  K  is  a  vector,  then  F  must  be  linear  in  all  the  vector  in 


the  list  (58),  plus  any  vectors  which  can  be  formed  in  an  invariant  way 

from  the  vectors  and  tensors  in  that  list.  Thus 

F  =  E  a.  V.  (59) 

1  1 

where  are  the  objective  vectors  so  formed.  Each  scalar  coefficient 

a^  can  be  a  function  of  all  the  scalars  and  scalar  invariants. 

If  F  is  a  second  order  tensor,  then 


F  =  E  B.  (60) 

where  T.  are  the  objective  tensors  which  can  be  formed  from  the  list 
J 

(58).  The  scalars  can  be  functions  of  all  the  scalars.  If  F  is 

symmetric,  we  need  only  consider  symmetric  if,. 
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The  general  approach  (Drew  and  Lahey  1979)  gives  a  problem  too  big 
to  be  practical  or  meaningful.  Thus,  we  shall  consider  here  forms  for 
the  stresses,  etc.,  which  we  expect  to  be  important  in  two-phase 
flows.  Obviously,  we  shall  have  the  constant  worry  that  something  has 
been  forgotten;  however,  the  alternative  of  stopping  at  this  point, 
awaiting  evaluation  of  over  two-thousand  scalar  coefficients  seems  less 
attractive. 

Stresses 

Let  us  now  discuss  specific  considerations  for  the  specific 
constitutive  equations  needed.  Let  us  begin  with  the  stresses.  The 
extra  stress  (viscous  or  elastic)  and  turbulent  stresses  are  combined  as 
+  o^.  Macroscopically,  at  least  in  the  context  of  this  model,  one 
sees  no  way  to  separate  them.  Microscopically,  one  is  due  to  the  actual 
stresses,  and  the  other  is  due  to  velocity  fluctuations  from  the  mean. 

In  many  cases  of  interest  the  extra  stress  is  smaller  than  the  turbulent 
stress,  one  reason  being  that  any  appreciable  slip  between  the  phases 
generates  velocity  fluctuations,  which  appear  in  this  model  in  the 
turbulent  stresses.  These  fluctuations  can  transport  momentum.  The 
conceptual  situation  is  analogous  to  that  for  single-phase  turbulent 
momentum  transport,  with  "eddies"  carrying  momentum  across  planes 
parallel  to  the  shearing  direction  and  losing  it  into  the  mean  flow. 

The  fluctuations  in  single-phase  turbulent  flow  are  there  because  of 
instabilities  in  the  laminar  shear  flow;  in  two-pha'^a  flow  they  also 
come  from  fluctuations  generated  due  to  the  flow  around  individual 


particles 
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Drew  and  Lahey  (1981)  have  generated  a  model  for  the  turbulent 
stresses  in  bubbly  air-water  flow.  The  model  they  deduce  in  this  case 
has  the  form 


O  =  2ui:D_  .  +  a. I  +  b  (v  -  V  )(v,  -  V  )  (61) 

2  22fU  2  21  2  1  2 

T 

for  the  liquid  phases.  Here  the  eddy  viscosity,  a2  is  an 

induced  eddy  "pressure",  and  b2  is  associated  with  bubble  passages. 

T 

The  parameters  y  ,  a  and  b,  depend  on  a,  the  bubble  radius  r, 

2  2  ^ 

the  relative  velocity  |v^  ”  ^2^  Euclidean  norm  of  D2  |j. 

Nigmatulin  (1979)  uses  a  cell  model,  assuming  inviscid  flow  around  a 
spherical  particle  in  each  cell,  to  calculate  values  of  a2  and  62* 

He  obtains 


b  =  —  —  p 
2  2  2 


(62) 


1  ,  ,2 
*2  =  6  *^2!^  -  ^2! 


(63) 


He  does  not  obtain  a  term  involving  the  eddy  viscosity  because  the  cell 
model  is  unsuited  to  obtaining  shear  flow  results*  His  results  are 
valid  for  dilute  flows.  Other  effects  are  present  in  most  turbulent 
flows. 


Drew  and  Lahey  (1981)  used  mixing  lengths  to  evaluate  the 
coefficients  in  eq.  (61).  The  data  available  thus  far  (Serizawa  1975) 
allows  only  evaluation  of  the  ratios  of  various  Karman  constants, 
nonetheless,  the  importance  of  the  problem  in  nuclear  reactor  technology 
has  spurred  the  acquisition  of  more  direct  data  in  special  cases  (Lance 
1981 ).  This  data  should  be  available  shortly,  and  should  allow  the 
evaluation  of  32  and  b2  for  bubbly  flows. 
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Evidence  (Serlzawa  1975)  also  indicates  that 


o 


1 


(64) 


and  is  negligible  in  dispersed  bubbly  flows.  A  model  for  analogous 

to  eg.  (61 )  is 

®1  =  2^1°1,b  *  *1^  ^‘^1  ■  ^2^^^  ■  ^2’  • 

If  the  particulate  phase  follows  the  fluid  phase  closely,  the  quantities 

in  equations  (61)  and  (65)  representing  both  velocity  scales  and  mixing 

lengths  of  large  eddy  processes  should  be  approximately  equal.  If  this 

is  the  case,  then 


M 


T 

1 


(66) 


The  remaining  terms  in  eg.  (65)  arise  at  least  in  part  due  to  the 
motion  of  the  particles  through  the  fluid.  No  constraints  are  placed  on 
these  coefficients. 

Ishii  (1975)  argues  that  for  the  viscous  fluid  phase,  the  average 
of  the  microscopic  viscous  stress  lead  to 


^2^0  2 


b( i-g) 
1-a 


[7o(i 


v^)  +  (. 


-  V2)Va] 


(67) 


where  is  ^be  exact  viscosity  of  the  fluid,  and  b  is  called  the 

mobility  of  phase  two.  He  argues  that  b(1-a)  ^  1  for  o  near  zero. 

No  experimental  evidence  has  been  offered  to  verify  (67). 

In  the  case  when  the  particulate  phase  consists  of  solid  particles, 
it  is  usually  assumed  that  3  0.  If  the  particulate  phase  is  also  a 
viscous  fluid,  then  it  is  sometimes  assumed  that 
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^  .  (68) 
This  corresponds  to  Ishii's  result  (see  eq.  67))  with  b  =  0. 

Murray  (1965a)  offers  a  model  for  a  particle-fluid  flow  in  which  he 
assumes  that  the  particles  experience  a  viscous  stress  with  a  viscosity 
coefficient  proportional  to  the  fluid  viscosity.  His  argument  indicates 
that  the  force  he  is  calculating  is  not  a  stress,  but  an  interfacial 
force. 

Pressure  Relations 

The  pressure  differences  i  discussed  next.  First,  we 

note  that  with  the  assumption  of  local  incompressibility  of  each  phase, 
we  must  constitute  all  but  one  of  p^,  P2,  p^  ^  and  P2  The  other 

will  then  be  a  solution  of  the  equations  of  motion. 

The  simplest  assumption  for  pressure  differences,  is  to  assume 
none.  That  is,  assume  Pjc  “  P)t  i  “  1,2.  This  assumes  that 

there  is  "instantaneous"  microscopic  pressure  equilibration  which  will 
be  the  case  if  the  speed  of  sound  in  each  phase  is  large  compared  with 
velocities  of  interest.  In  applications  which  do  not  deal  with  acoustic 
effects  or  bubble  expansion/contraction,  this  assumption  is  adequate. 

For  a  discussion  of  bubble  oscillation  effects,  see  the  article  by  van 
Wijngaarden  (1972).  Cheng  (1982)  studies  the  consequences  of  several 
modeling  assumptions  on  wave  propagation  in  a  bubble  mixture. 

In  situations  where  surface  tension  is  important  and  no  contacts 
occur  between  the  particles,  we  assume  that  the  jump  condition  (  4) 
becomes 

1  2 

p  -  p  = 


OX 


(69) 
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where  <’  is  the  exact  curvature  of  the  interface.  Multiplying  by 


and  averaging  gives 

(P,  •  -  P,  • =  ax  Vo  ,  (70) 
Where  <  is  the  average  mean  curvature*  Since  this  equation  holds  for 
all  Va,  we  have 


p  -  p  =  ax 
*^1,1  2,i 


(71  ) 


If  the  surface  tension  is  negligible,  equation  (2.72)  gives  Pi^i  = 


P2,i‘ 

When  contacts  occur  between  particles,  these  considerations  are  not 
always  simple.  The  contacts  provide  mechanisms  for  causing  the  average 
pressure  at  the  interface  in  one  phase  to  be  higher  than  the  other. 
Consider  solid  particles  submerged  in  fluid,  with  no  motion,  and  with 
the  particles  under  a  compressive  force  (supporting  their  own  weight,  or 
the  weight  of  an  object,  for  example).  The  resulting  stress  will  l>e 
transmitted  from  particle  to  particle  through  the  areas  of  particle- 
particle  contact.  Just  on  the  particle  side  of  the  interface  where  the 
contact  occurs,  the  stress  (pressure)  level  will  be  large.  At  places 
where  no  contact  occurs,  the  stress  will  be  equal  to  the  pressure  in  the 
surrounding  fluid.  Under  normal  circumstances,  the  contact  areas  will 
be  a  small  fraction  of  the  total  interfacial  area.  Thus,  the 
approximation  i  ^  ^2  i  ”  ^i  usually  valid.  In  the  more  general 
situation,  these  contacts  may  be  intermittent,  when  the  particles  are 
bumping  together,  or  may  be  constant,  as  when  the  particles  are  at 
rest.  In  this  case,  also,  we  expect  i  “  ^2  i  ~  Pi* 
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Stuhmiller  (iy77)  argues  that 


P2  -  P2,i  =  ^‘^2l''2  - 


(72) 


where  C  is  taken  to  be  a  constant.  The  argument  leading  to  (72)  is 
exactly  the  calculation  of  form  drag,  which  is  usually  part  of  the 
interfacial  drag.  It  is  not  clear  whether  using  (72)  and  a  drag  law 
(79)  includes  the  form  drag  consistently. 

We  shall  assume  that  for  the  fluid  phase,  P2  =  Pj^«  To  allow  for 
the  possibility  of  the  extra  stress  due  to  contacts,  in  the  particulate 
phase  when  it  consists  of  solid  particles,  we  shall  write 


Pi  =  Pl,i  Pc  =  P2  Pc 


(73) 


where  p^  is  the  pressure  in  the  particles  due  to  contacts.  A  model  is 

needed  for  p^.  Several  modelers  (Gough  1980,  Kuo  et  al.  1976)  have 

used  a  model  with  p_  =  p  (o),  with  p  (o)  =  0  for  a  <  o  ,  where 

c  c  c 

01  is  the  packing  concentration  of  the  particles,  and  p  (o)  rapidly 
c  c 

increasing  for  a  >  This  further  suggests  an  "incompressible" 

model,  with  p^  =  0  for  ot  <  a^,  and  a  =  a^,  and  p^  unspecified, 
otherwise.  In  this  latter  region,  the  particle  momentum  equation  is 
meaningless.  Such  a  model  has  not  been  exploited. 

_Interfacial  Momentum  Transfer 

Let  us  now  consider  the  interfacial  momentum  transfer  terms.  In  a 

sense,  the  crucial  modeling  for  two-phase  flow  is  done  with  these  terms, 

since  they  are  the  terms  which  couple  the  two  phases  together. 

Let  us  consider  M^.  In  analogy  with  Ishii  (1975),  we  write 

H  =  a<  Va  (74) 

m 

where  <  is  the  average  mean  curvature  of  the  interface.  In  this  case. 
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the  momentum  jump  condition  becomes 

+  p  ^  -  p  Va  =  OK  7a  .  (75) 

Using  condition  (71),  we  have 

=  0  .  (76) 

Therefore,  we  need  only  provide  a  model  for  then  is  determined 

by  the  jump  condition  (76). 

The  interfacial  momentum  transfer  contains  the  forces  on  the 

particulate  phase  due  to  viscous  drag,  waJce  and  boundary  layer 
formation,  and  unbalanced  pressure  distributions  leading  to  lift  or 
virtual  mass  effects,  except  for  the  mean  interfacial  pressure.  The 
model  which  we  discuss  here  should  contain,  as  much  as  possible,  the 
above  forces.  Indeed,  this  motivates  the  inclusion  of  certain 
quantities  in  the  list  of  variables  (58). 

We  postulate 


+  A,(v_  -  V, )*D  +  A  (v  -  V,  )*D  +  A  (▼  -  v  )*7(v  -  v  ) 

3  2  ’  I  /  D  42  I  2^D  ^2  1  2  1 


where  A^  -  Ag  are  scalar  functions  of  the  invariants. 


(77) 


The  first  term  represents  the  classical  drag  forces.  To  conform 
with  customary  useage,  we  write 


^  =1“  "^2  7-  -  ^2! 

where  r  is  the  effective  particle  radius,  and  is  the  drag 

coefficient.  It  is  usually  assumed  that  =  Cp(a,Re),  where 
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Re  =  ~  ^2  particle  Reynolds  number.  A  careful 

study,  including  extensive  comparisons  with  data,  is  given  by  Ishii  and 
Zuber  (1979).  Their  conclusions  for  the  drag  coefficient  are  summarized 
in  Table  1 ,  which  is  adapted  from  their  paper. 

The  combination 


3v  3v 

is  an  acceleration  term.  Drew,  Cheng  and  Lahey  (1979)  write 

""2  =  “  ' 

^5  "  A2(1-X(a))  , 


(79) 

(80) 
(81  ) 


1 


where  is  referred  to  as  the  virtual  volume.  If  and  no 

spatial  velocity  gradients  are  present,  then  (79)  reduces  to  the 
classically  accepted  virtual  mass  force.  It  is  less  easy  to  obtain 

correlations  for  than  for  C^,  since  the  ratio  of  virtual  mass 

* 

forces  to  drag  forces  is  of  order  V  t  /r,  where  V  is  a  velocity 
scale  and  t  is  a  time  scale.  In  order  for  the  virtual  mass  force  to 
be  appreciable,  the  time  scale  must  be  of  the  order  r/V.  For  most 
multiphase  flow  applications,  r  is  small,  and  V  is  not.  Thus,  we 
see  virtual  mass  effects  only  at  relatively  high  frequencies. 

Zuber  (1964)  considered  finite  concentration  effects  on  virtual 
mass  by  considering  a  sphere  moving  inside  a  spherical  "cell".  He 
obtained 

1  1+7«  1  1 

(82) 


^  .  1  1+20  1  3 

c  (o)  - ^  a 

VM  =  2  1-0  =2  2 


for  small  o.  Nigmatulin  (1976)  also  calculates  the  virtual  mass 

coefficient  from  a  cell  model.  He  obtained  C  (o)  ^  —  to  0(o).  The 

VM  =  2 

form  of  the  acceleration  which  he  derives  is  not  objective.  It  is 
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possible  that  this  is  due  to  the  inability  of  cell  models  to  adequately 

deal  with  velocity  gradient  terms.  The  value  of  the  virtual  mass 

coefficient  derived  by  van  Wijngaarden  (1976)  is 

=  0-5  +  1.39a  .  (83) 

VM 

Moheyev  (1977)  obtained  an  empirical  correlation 

C  „  (a)  =  0.5  +  2.  1o  (84) 

VM 

by  using  an  electrodynamic  analog.  Thus,  there  seems  to  be  a  concensus 

that  C  -*•  "  as  a  +■  0.  There  is  less  agreement  on  the  0(a) 

VM  2 

correction.  In  addition,  there  are  several  models  in  the  literature 
(Wallis  1969,  Hinze  1972)  which  use  non-objective  forms  of  the  virtual 
mass.  There  is  no  compelling  evidence  at  this  point  as  to  whether  the 
objective  or  non-objective  forms  conform  more  closely  to  reality. 

Cheng  (1982)  has  done  extensive  predictions  using  virtual  mass 
models  for  accelerating  air-water  bubbly  flows.  For  low  frequency  small 
disturbances  in  a  one-dimensional  flow,  he  found  that  the  value  of 
Cym  had  a  significant  effect  on  the  propagation  speed,  but  that  the 
data  were  scattered  so  that  a  particular  value  could  not  be  estimated. 
For  higher  frequency  waves  (sound  waves)  and  for  critical  flow, 
compressibility  effects  were  important.  For  nozzle  flow,  the  virtual 
mass  effects  were  insignificant.  None  of  the  calculations  showed  any 
effect  of  different  values  of  X. 

The  term 

\^^2  -  ^1>'''2,b 

contains  the  effect  of  the  lift.  We  write  A,  =  a  L.  On  the  other 

4 


hand,  the  forces  represented  by 
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A  (v,  -  v_)*[VCv  -  V  )1^  have  no  analogs  in  single  particle 
6  12  12 

calculations,  and,  no  observations  confirm  their  presence.  Thus  we 
assume  ^2  ~  ~  Very  little  experimental  evidence  for  exists. 

Interfacial  Geo_metr_y  Models 

It  is  evident  from  eqs.  (71)  and  (78)  that  a  relation  is  needed  for 

expressing  the  geometry  of  the  interface.  Several  such  relations  have 

been  used;  however,  none  is  so  compelling  as  to  be  called  general.  The 

one  which  comes  closest  is  due  to  Ishii  (1975)  (derived  for  the  time 

average;  the  generalization  is  obvious).  We  shall  discuss  this  later. 

First,  we  present  the  simplest  models.  If  the  particles  are  solid 

spheres  each  of  the  same  radius  r,  we  have 

o  =  Y  n  r^  n  ,  (86) 

where  n  is  the  number  density.  Other  quantities  of  interest  are 

immediately  obtainable  from  a  and  r,  for  example,  the  average 

3a 

inter facial  area  per  unit  volume  is  Note  that  n  is  conserved 

since  a  is;  specifically,  we  have 

j-  +  V.n  =  0  .  (87) 

If  the  particles  are  not  mon'idisperse,  it  may  be  necessary  to  derive  a 
model  where  each  size  of  particle  is  treated  separately.  Some  recent 
efforts  along  these  lines  show  great  promise  (Greenspan  1982).  Particle 
breaking,  agglomeration,  or  accretion/erosion  also  necessitates  a  more 
complicated  model.  An  effort  (without  macroscopic  mechanical 
considerations)  is  given  by  Seinfeld  (1980).  Both  approaches  are 
essential  geared  to  an  equation  analogous  to  (37)  for  a  number  density. 
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A  more  general  approach  is  suggested  by  Ishii  (1975)  (for  time 

averaged  equations),  where  he  writes 

ds 

+  y^V^s  =  ,  (88) 

where  s  is  the  interfacial  area  per  unit  volume,  which  in  the  present 
notation  is  <3x/3n>,  is  the  averaged  interfacial  velocity,  and 
is  a  source  of  interfacial  area,  due  to  "bulging”  of  the  interface. 

This  approach  is  probably  loest  suited  to  flows  where  the  particulate 
phase  (bubbles  or  drops)  can  coalesce  or  brealc  up.  It  has  not  been 
investigated  to  any  degree. 
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model  CONSIDERATIONS 


The  main  aspect  of  two-phase  flows  which  makes  them  fundamentally 

different  from  a  single-phase  flow  is  the  interfacial  interaction  terms, 

and,  thus,  in  the  absence  of  interfacial  mass  transfer  the  difference  is 

d 

contained  in  the  interfacial  force  The  modeling  of  this  term  must 

proceed  cfirefully,  with  attention  paid  to  any  possible  inconsistencies 
or  untoward  predictions.  A  canonical  problem  involving  simplified  two- 
phase  flow  models  has  surfaced  in  the  literature;  namely,  that  the  one¬ 
dimensional,  incompressible,  inviscid  flow  equations  without  virtual 
mass  effects  are  ill-posed  as  an  initial  value  problem  (Ramshaw  &  Trapp 
1978). 


To  see  this,  consider  the  equations  with  C 


VM 


0  and  L  =  0. 


(This  represents  the  model  used  in  most  practical  calculations.  ) 


3a  „ 

3t  'Tx"  "  ° 


(89) 


3a  , 

Ft  ■  3x""  -  ° 


(90) 


3 


3v. 


3x“^  =  Tx-  "  ^'"2  -  ''l^ 


(91  ) 


9  V  3v 

(1-a)P2(3--^-  +  v^  y/)  =  -d-a)  +  A,(v^  -  v^) 


(92) 


If  the  pressure  is  eliminated  from  system  (89)  -  (92)  and  the 


3u 


3u 


resulting  equations  are  expressed  in  the  form  A  +  B  =  C,  where 


u  =  [a.v^.v^l  ,  A  and  B  are  3  3  matrices  and  C  is  a  3-vector, 

dz . 

then  the  characteristics  of  the  system  are  given  by  ,  i  =  1,2,3, 


I 
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where  v,  are  the  three  characteristic  values  of  det[Av-B]  =  0. 

1 

Note  that  since  C  contains  no  derivative  terms,  it  is  not  involved  in 
determining  the  characteristics.  The  characteristic  values  for  the 

system  as  given  are  a  pair  of  complex  conjugate  roots  v  _  =  ot  ±  iB 

\  t  ^ 

(with  B  0),  and  =  “  {resulting  from  the  assumption  of 
incompressibility).  This  implies  (Garabedian  1964)  that  the  initial 
value  problem  is  ill-posed. 

One  implication  of  the  ill-posed  nature  of  the  system  can  be  seen 
as  follows.  Consider  the  linear  stability  of  a  constant  solution.  The 
solution  of  the  linearized  equations  has  the  form  u(x,t)  = 

Ug  +  u^e  ,  where  k  is  real.  Substitution  into  the  linearized 

system  results  in  the  eigenvalue  problem 

[A(Uq)\)  +  B(uQ)ik  -  C(Uq)]u^  =  0  (93) 

for  V.  The  eigenvalue  equation  is 

iv 

det[A(uQ)  - - k — 

For  k  large,  two  roots  of  (94)  are 

=  a  ±  iB  .  (95) 

Thus, 


v,  -  *  k(±B  -  ai)  .  (96) 

As  long  as  B  0,  one  of  the  w's  will  have  positive  real  part, 
indicating  linear  instability.  Note,  too,  that  as  k  increases,  the 
growth  rate  |Bk|  increases.  This  shows  that  finer  disturbances  grow 
faster.  This  implies  that  at  any  t  >  0,  maxju^e  |  can  be  made  as 

large  as  desired  by  taking  k  sufficiently  small.  Contrast  this  with 
the  behavior  of  an  instability  of  a  well-posed  system,  where  the 
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^  fC  X 

solution  u^e  can  be  made  as  small  as  desired  at  a  finite  t  by 

making  u^  small,  for  any  k.  The  instability  in  the  well  posed  system 
has  some  realistic  physical  meaning,  while  the  instability  always 
present  in  the  ill-posed  system  suggests  that  the  model  is  not  treating 
small  scale  phenomena  correctly. 

The  extension  of  the  argument  to  three  space  dimensions  is 
straightforward,  and  shows  that  the  system  is  ill-posed  in  three 
dimensions  also. 

The  ill-posed  nature  of  the  model  makes  it  desirable  to  seek 
physical  modifications  of  that  simple  model  to  find  a  well  posed 
model.  An  effort  has  been  made  to  determine  whether  and/or  to  what 
extent  virtual  mass  makes  the  inviscid  model  well  posed.  Lahey,  et  al. 
(1980)  have  shown  that  the  model  becomes  well-posed  for  certain  values 
of  Cynj.  The  value  of  Cyj^  which  makes  the  system  well-posed  is 
usually  large  unless  a  is  extremely  small,  or  other  somewhat  non¬ 
physical  assumptions  are  made.  Thus,  even  though  it  is  a  possible  "fix" 
to  the  question  of  the  well  _vs^  ill-posed  nature  of  the  model,  virtual 
mass  does  not  seem  to  be  the  total  answer. 

We  note  here  that  if  the  viscous  and  Reynolds  stresses  are  included 
in  the  equation,  and  they  are  assumed  to  depend  on  the  equations 

become  parabolic,  and  hence  are  well-posed  if  enough  boundary  conditions 
are  used.  Also  Prosperetti  and  van  Wijngaarden  (1976)  have  shown  that 
using  certain  compressibility  assumptions  in  gas-liquid  systems  gives 


real  characteristics.  Stuhmiller  (1977)  shows  that  his  model  including 
the  form  drag  (see  eq.  72)  gives  real  characteristics  for  sufficiently 
large 
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The  question  of  the  role  of  the  inviscid  model  remains.  In  single¬ 
phase  fluid  mechanics,  the  inviscid  model  is  extremely  important,  and 
governs  the  flow  outside  of  boundary  layers,  shear  layers  and  other 
singular  structures.  In  two-phase  flows,  the  inviscid  equations  are 
ill-posed,  possessing  the  unrealistic  feature  that  small  scale  phenomena 
grow  rapidly.  This  seems  to  imply  that  an  inviscid  model  is  nowhere 
valid  for  two-phase  flows,  that  is,  viscous  or  eddy  stresses  are 
important  everywhere.  An  important  unanswered  question  is  whether  the 
inviscid  limit  of  the  two-phase  flow  equations  is  meaningful. 

Let  us  examine  another  possible  hypothesis  about  the  interfacial 
force  which  leads  to  an  interesting  form  of  the  virtual  mass  and  lift 
terms.  The  hypothesis  is  that  the  interfacial  force  must  be  dissipative 
for  all  conceivable  flows.  This  requirement  can  be  derived  from 
considerations  involving  an  entropy  inequality.  The  argument  is 
complicated  (Drew  and  Lahey  1981),  and  will  not  be  presented  here.  It 
depends  on  the  correct  grouping  of  several  terras  in  the  energy  and 
entropy  equations. 

If  we  form  the  total  kinetic  energy  equation  by  dotting  (42)  with 
1^,  and  adding  for  k  =  1,2,  we  obtain 
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2  2 


k=1 


k=l 


(97) 


The  terms  on  the  left  represent  the  rate  of  change  total  (mean 

flow)  kinetic  energy.  The  pressure  terms  contribute  to  the  stored 

2 

internal  energy.  The  term  1  **k^\^**k^*\  working  of  the 

t  d 

stresses.  The  two  terms  -  Z  ("r, +®,  ):^v,  and  -M  ‘(v.-v^)  are  the 

k  k  k  k  112 

(bulk)  dissipation  and  the  interfacial  dissipation,  respectively.  We 

require  the  total  dissipation  to  be  expressible  as 

3a 

+  V*B  +  C  ,  (98) 

where  A  >  0  and  C  >  0  for  all  possible  flows.  The  term  A 
represents  the  extra  kinetic  energy  attributable  to  interactions,  B 
represents  the  extra  kinetic  energy  flux,  and  C  is  the  (actual)  total 
dissipation. 

We  shall  present  the  results  of  the  disspativity  requirement  for 
=  0  and  =  0.  If  we  then  take  eq  (78)  for  M^,  and  assume 
that  Cynj,  X  and  L  are  constants,  we  have 


3 

^  3I 


t  "  (1-X)(v^-v^)]— ^  3^---} 


(V1-V2) 

- 2 - f  °^VMP2  ( ^-2  )  ( V,  -v^  )  ] 
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(99) 


If  we  examine  (99)  for  uniform  fields,  we  obtain  ttie  first  result  of  tlie 
dissipativity  requirement,  namely  that  >  0.  Furthermore,  we  note 
that  since  the  other  terms  depend  on  spatial  gradients,  and  we  can 
impose  them  at  will,  we  have 

X  =  2  ,  (100) 


^  =  -V2  • 

Note  also  that  the  requirement  that  A  >  0  results  in  >  0. 

d 

resulting  form  for  is 


(101) 

The 


(102) 


In  this  form,  the  virtual  mass  terms  are  in  conservation  form,  with  the 


convective  velocity  v^.  The  non-objective  nature  of  these  terms  must 
be  compensated  for  by  the  lift  term,  resulting  in  the  last  term  in 
(102).  Note  that  if  **2  ”  ^  ^  ^^2  fluid  vorticity,  thelast  term 

becomes  ®C  (v  -▼  )  *  «  .  Thus,  the  remaining  part  of  the  lift  is 
always  perpendicular  to  the  slip  and  the  fluid  vorticity. 
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